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ABSTRACT 


Analytic solutions for non-linear differential equations are 
difficult, time consuming and largely impractical for reasonably large 
values of the independent variable. The purpose of this study was to 
develop a technique for analytic (algebraic) solutions of autonomous and 
nonautonomous equations of the type 

q n 

y ] C n - ^1T = + kN etc “] 

n=0 dt 


with the help of a digital computer. 

In the equation n is an integer and q is a small integer. c n is 

a constant. g(t) is a forcing function. N jf, f', f " etclj is the non-linear 

term while k is the usual "small" parameter. N does not contain the 

independent variable t (time) explicitly. f(t) is a continuous bounded 

function with finite initial conditions. 

Two operational transform techniques have been programmed for 

the solution of equations of this type. To develop computer techniques 

only relatively simple non-linear differential equations have been con- 

/ 

sidered. The theorem of Poincare assures the convergence in all cases 
considered for "small" values of k. 

In the few cases considered-it has been possible to assimilate 
the secular terms into the solutions. 

For cases where f(t) is not a bounded function a direct series 
solution is developed which can be shown to be an analytic function. 

All solutions have been checked against results obtained by 
numerical integration for given initial conditions and constants. 


While the results of this study may be regarded as experimental 
in character it seems evident that at least certain types of non-linear 
differential equations not only can be solved with the help of a digital 
computer but that except for quite elementary equations must be solved 
in this way. 
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1. INTRODUCTION 


Engineers and physicists have been interested in non-linear 
problems for many years. While graphical methods or hand-calculated 
values for numerical integration could sometimes be employed, the use- 
fulness of such solutions was largely restricted by the tedious calculations 
and techniques. With the advent of digital and analog computers numeri- 
cal and graphical solutions could be obtained for most engineering prob- 
lems. Such solutions require definite numerical values for design con- __ 
stants and initial conditions and many different runs are necessary to 
establish the character 'of the solution. 

In many cases analytic (algebraic) solutions might prove quite 
helpful if not too complicated since the effect of changing design con- 
stants or initial conditions could be recognized from the form of the 
solution. This would be particularly true if the results could be approxi- 
mated with two or three terms. Unfortunately analytic solutions for non- 
linear differential equations are both difficult and time consuming and, 
in most cases, quite impractical for large values of the independent" vari- 
able. 

During the past seven years the author and his students conducted 

a series of investigations leading to analytic solutions for non-linear and 

linear differential equations using a digital computer. Madhu 1 showed 

that it was possible to solve a non-linear differential equation using a 

power series in the independent variable, time. This approach was both 

cumbersome and unsuited to machine manipulation. It was evident that 

if machine calculations were to be used an operational or transform type 

of solution would be required. 

2 

Crecraft developed solutions for several non-linear equations 

3 

but only for cases where the solutions were known. Shiva showed that 

a linear differential equation with a variable coefficient could also be 

solved by a technique which could be adapted to computer operations. 

4 

Regan programmed a computer for certain derivative operations and 
5 7 8 

Erk'anli , Cheng and Wang developed tables of transforms for special 
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functions. Numerical solutions for the van der Pol equation were ob- 
0 

tained by Chen while numerical solutions for the Mathieu equation were 
9 

computed by Lee „ 

The purpose of this report is to present the results of an experi- 
mental study in the use of a digital computer to obtain analytic solutions 
for certain types of non-linear differential equations. The solutions are 
limited to bounded analytic functions and convergence is assured by 
Poincare! Theorem . Appropriate cases include many of the equations 
usually discussed in text-books such as certain forms of the van der Pol 
and Duffing equations. 

One section is devoted to direct series solutions where the results 
are unbounded analytic functions. 

It is sometimes said that calculations which can be done by hand 
can always be programmed on a digital computer. Unfortunately it is 
not unusual to find that the whole procedure must be revised to carry 
through a given problem. In the case of non-linear solutions the power- 
ful perturbation method would be very difficult to program due to the 
series of solutions required for the unknown coefficients. An operational 
type of solution seemed advisable since a computer could easily be pro- 
grammed to perform the required operations. It should be emphasized 
that the Laplace Transform is employed only in evaluating the equation 
as presented. In particular the solution of several simultaneous non- 
linear differential equations by Laplace Transforms is not involved. 

Two different techniques are described. The first type is closely 
associated with the development of a Maclaurln 1 s series by repeated 
differentiations while the second type evaluates the time function directly, 

A discussion of the conditions for convergence of the operational 
solution will be followed by an example presented in some detail. 

An outline of the computer operations is then followed by a number 
of examples which show the power and weakness of the method. In each 
case the numerical values are substituted in the analytic solution to obtain 
a graphical solution which is then compared with a solution obtained by 
numerical integration. 
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2. TYPE OF EQUATION CONSIDERED 


The primary purpose of this study was to develop a method for 
obtaining analytic solutions for non-linear differential equations using 
a digital computer. For this reason the equations selected were as 
simple as possible. 

All equations are of the form; 


q 


c 

n=0 


,n 


n dt n 


f(t) = g(t) + kN f, f’, f"', etc.J 


( 1 ) 


Here n is an integer and q is a small integer. c n is a constant. g(t) is 
a forcing function. N jf, f', f" , etc.] is the non-linear term in polyno- 
mial form while k is the usual n small ,r parameter. The equation is 
autonomous if the forcing function is a constant and nonautonomous if it 
is not a constant. f(t) is a bounded continuous single valued analytic 
function. All terms are limited to functions which may be developed in 
Taylor's series. 

Equation (1) is the form usually considered in text books on non- 
linear differential equations and convergence of the solution is assured 
if the parameter k is small enough on the basis of the theorem of 
Poincar^^. 

All solutions are therefore analytic functions. 



3. CONVERGENCE OF THE OPERATIONAL 
SOLUTION 


4 


Consider a function of the type 
00 

f(t) = > s k f (t) (2) 

/ . mm ' 

m=0 


f(t) is a bounded, continuous, single valued function uniformly convergent 
in some closed interval 0 < t < h, h > 0. All f (t) are similarly defined 
(See Reference 11). In addition 



0 


a restriction required by the technique in the solution. 
The Laplace direct transform is 


F(p) = p 



f(t) e pi: dt 


(3) 


where p is a complex number with real part greater than zero. 

Equation (3) is to be carried through by integration by parts. 
Since all f m (t) are bounded functions the upper -limit contributes nothing 
to the final result: 


E 

m 


(p) = P 



f (t) e" pt dt = 
m 


oo 


n=0 



( 0 ) 


n 

P 


(4) 


The result in p n is always absolutely and uniformly convergent in p 
(See Reference 12). 

If the upper limit on the integral in (4) is removed the '.result is an 
"open-end integral" and the result in (4) may then be obtained by a series 
of differentiations carried through by a computer. 
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When all f (t) converge uniformly and the integral in (4) con- 


verges 


■/' 


f (t) e pt dt = p 
m v 


of 


00 


f n (0) -i-j- e~ pt dt 
m n l 


23 

n=0 


oo f n (Q) 
m 


n 


= F (p) . 
. m ^ 


( 5 ) 


That is, the function and the series in t n have identical transforms and 
all F m (p) are absolutely and uniformly convergent; f (t) is in terms 
of a Maclaurin's series. 

Now uniformly converging series can be summed (Reference 11): 


00 


F(p) = 



00 k f n (0) 
m m 



n 


( 6 ) 


m=0 n=0 p 

If (6) is now inserted in the inverse transform 

C+joo 


f(t) = 


/ ' C+Jo 
C“jco 


e P‘dp 


( 7 ) 


-n. 


each of the p th terms may be evaluated separately and the results added 
to obtain 


oo 


f(t) = 


23 

m=0 n=0 



°° k f n (0)t n 
mm 


n i 


( 8 ) 


Data' comes from ’the- machines innthe form 


F.'(p)/%/ 



k f n (0) 
mm 

n 


( 9 ) 


m= 0 n= 0 

where q and r are integers. Normally q < 6 and r < 30. The symbol A<# 
indicates an approximation. 
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The direct series result equivalent to (8) is then 


q r 

«*> ~T2 

m= 0 n= 0 



k f n (0)t 
mm 

n I 



( 10 ) 


The series (8) is a convergent series by definition and since all 
f m (t) are uniformly convergent (10) is always the first portion of a con- 
verging series which may be carried as far as necessary. 

If the raw data of (9) are sortedon k and the p n th terms summed 

m 

in terms of known functions, 

q 

F(p)/^ j: 

m=0 

If these F (p) are now evaluated by the inverse transform, 

q 

< 12 > 

m=0 

Since f(t) and all f (t) are bounded functions and the series is 
uniformly convergent, f(t) converges if q = oo. 

Equation (12) is a "summation solution". 

Such summation solutions are much more effective than the direct 
series solutions in producing results. 

In the direct series solution all data are evaluated directly. It is 
not necessary to recognize the F m (p) in order to obtain a direct series 
solution. 

On the other hand, a summation solution uses a few terms to 
recognize an operational term and the resulting solution converges to a 
usable result for much larger values of t than the direct series type of 
solution. 

As might be expected summation solutions quite frequently have 
secular terms. These terms can be assimilated 'into the solution in some 
cases. With a computer presumably a sufficient number of them could be 
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obtained so that the series thus formed would be a good approximation 
for the unknown function over the desired range. In this case the con- 
dition of uniform convergence is not usually valid. 

All summation solutions do not have secular terms, and in at 
least one instance, the solution for the pendulum, the frequency of oscilla- 
tion can be predicted and no secular terms appear. This suggests that 

with proper summation of terms A (f{0), f'(0) etc., k 

} .n 

where the coefficients A and the time functions contain k and 

n 

the boundary conditions the results are the same as for Lindstedt's 
method^. 

Operational solutions for equations defined in Section 2 are formal 
solutions of the equations. Convergence is assured for bounded functions 
if the parameter k is small enough and a sufficient number of terms is 
available . 

In the special case where the formal solution f(t) can be recognized 
as a bounded, continuous, single valued function uniformly convergent in 
some closed interval 0<t<h, h>0 and k is small enough the result con- 
verges. Such results are analytic functions. 
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4. DERIVATIVE OPERATIONAL METHOD 


Two operational methods have been developed for the solution of 
suitable non-linear differential equations using a computer. The 
"derivative" or "open-end integral" method is explained in this section 
by the solution of a well known first order equation. Section 5 explains 
the operations for the " expansion" method. 

Consider the equation 

f 1 + Af = -Bf 3 (13) 


where A and B are positive constants. B is "small" . Physically one. 

* 

may visualize this equation as non-linear braking for a rotating mass. 
The solution f(t) is therefore a bounded function. 

Rewrite (13) as 


f = -Af - Bf 3 

and let f(0) = a, the only initial condition. 
The direct transform is 


F(p) = p 



f(t) e" pt dt. 


(14) 


(15) 


Integrating by parts: 


F(p) = p 


-fe' pt 

P 

oo 

>- 1 1 & 
+ 

/” 00 

f' e“ pt dt 


0 

■J 

0 



V 



( 16 ) 


Since the function f(t) is bounded the upper limit in (15) never 
contributes to the final result and the integration by parts turns into 
a series of differentiations with the lower limit f(0) = a inserted. 

Continuing the integration by parts after substituting (14) in (16) 
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F(p) = a + 



(-Af - Bf 3 ) e' pt dt 


s CO 


= a - A j 


0 


f e 





= a 


Aa ^ _A_ 
P P 


r 


oo 


f e' pt dt 


Ba" 


J 


r 


00 


3f 2 f» e” pt dt. 


(17) 


J 


0 


This procedure can easily be programmed .on a computer using 
a defined operation which may be designated as the "open-end integral": 


f(t) e" pt dt 




n=0 


f n (0) 

n 


(18) 


JO. 


where f (t) is the n th derivative. 

If the inverse transform is used term by term 
oo 

'X ^ f^O) t n 

= x , 

n=0 


( 19 ) 


which is Maclaurin's series for f(t). 

Continuing with the derivative operations with Equation (18) the 
computer gives, for the first seven derivative operations: 
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X aA + a 3 B . aA 2 + 4a 3 AB + 3a 5 B 2 

rip; - a - — + 9 


aA 3 + 13a 3 A 2 B + 2 7a 5 AB 2 + 15a 7 B 3 


A Q Q ^99 7 S Q 4 

aA + 40a A B 4- 174a A B 4- 240a AB + 105a B 


r aA + 121a A B + 990a A B + 2550a A B + 2625sT AB 
L+ 945a U B 5 


aA 6 + 364a 3 A 5 B + 5313a 5 A 4 B 2 + 22,880a 7 A 3 B 3 
+ 41,475a 9 A 2 B 4 + 34,020a 1;l AB 5 +- 10,395a 13 B 6 


aA 7 + 1093a 3 A 6 B + 27,657a 5 A 5 B 2 + 186,165a 7 A 4 B 3 
+ 532,875a 9 A 3 B 4 + 747,495a U A 2 B 5 + 509,355a 13 AB 6 
_+ 135,135a 15 B 7 


+ + + + 

The Maclaurin's series equivalent to (19) is: 


( 20 ) 


(aA + a 3 B)t , (aA 2 + 4a 3 AB + 3a 5 B 2 )t 2 
+ yr 

(aA 3 + 13a 3 A 2 B + 27a 5 AB 2 + 15a 7 B 3 )t 3 
3 : 


f(t) = a - 
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, (aA 4 + 40a 3 A 3 B + 174a 6 A 2 B 2 + 240a 7 AB 3 + 105a 9 B 4 )t 4 

+ 4 1 

“3 A k q o 79*1 q A 

- (aA + 121a A B + 990a A B + 2550a A B + 2625a AB 
+ 945a i;L B 5 )t 5 

5 l 

+ (aA 6 + 364a 3 A 5 B + 531a 5 A 4 B 2 + 22 a 800a 7 A 3 B 3 
+ 41,475a 9 A 2 B 4 + 34,020a 11 AB 5 + 10 J 395a 13 B 6 )t 6 

- (aA 7 + 1093a 3 A 6 B + 27,657a 5 A 5 B 2 + 186, 165a 7 A 4 B 3 

+ 532,875a 9 A 3 B 4 + 747,495a U A 2 B 5 + 509,355a 13 AB 6 
1^77 

+ 135,135a °B )t 

7 

+ + + + ( 21 ) 

This is the direct series solution for Equation (13) approximated 
by the first eight terms. The (infinite) series converges if B is "small 
enough". (However, see Section 11.) 

Figure (1) shows, in graphical form, the result of substituting 
the indicated values in nineteen terms of Equation (21) extended. The 
values for the direct series solution are shown superimposed on a 
solution by numerical integration. 

Finite power series solutions of this type only represent the 
functions for relatively small values of the independent variable. Eventually 
the sign of the last term causes the series solution to go above or below the 
numerical solution. 




3 * + 1.0 

B- + O.I 
A“ + 0.2 

Numerical Evaluation 

xxx Direct Series Solution 



to 
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If the computer data for Equation (20) (nineteen terms) are sorted 
on B n and the p n terms are summed (if } pjj is large) the first three 
terms are 


F(p) = 


a Bp 


3a 5 B 2 p 


(p+A) ‘ (p+A) (p+3A) (p+A) (p+3A) (p+5A) 


+ + + • 


If these terms are evaluated by the inverse transform 
3 

-,,, -At a B -At ,, -2 At, 

f(t) = a e sta — e (1 - e ) 


2A 

3 a 5 B 2 -At 

-r e 


1 1 -2 At, 

(1 - e ) 


( 22 ) 


5 a 7 B 3 -At „ -2 At,' 

16 -p— e (1 ’ e > 


4 - + + + . 


This is the "summation solution". , 

Convergence is assured if B is small enough. 

To test for uniform convergence by the " M" test 


(23) 


CO 


12 M n = 


„ , a S B ,3 

- a + -23~ f T 


5^2 , 73 

s. B 5 sl B 

“72 + 16 “73“ 


n=0 


(24) 


since all of the exponential coefficients of (24) are less than one for any 
value of t since 


00 



M = 


n=0 


n il 


a 2 B 

A 


(25) 


provided B is small enough. 

The series (23) is therefore absolutely and uniformly convergent 


t > 0 for 


0 < < 1, 


(26) 
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In this case the terms of Equation (23) may be summed by binomial 
expansion provided 

(1 - e" 2At ) < 1 (27) 

giving 



If the solution of one case is known the result may sometimes be 
stated in a more general form. 

The solution of the equation 

- + Af + Bf q = 0 (29) 


where q is a positive integer is: 



for suitable values of the constants. 


( 30 ) 
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5. SOLUTION BY EXPANSION 


The Laplace Transform method of solution may easily be adapted 
to machine calculations using a method which has been designated as the 
"expansion” method since the key operation is to expand the solution as 

f = f x + f 2 + f 3 + + + + . (31) 


The results appear to match solutions obtained by the summation 
derivative process of Section 4. In fact the operational root functions 
for the summation procedure are readily predicted by a rough study 
using the expansion principle. 

As an example in the use of the method consider the equation 


f ' + xf = ge ^ + k (. f 2 -' f ff * ) 


(32) 


with f(0) = a. f(t) may only be a bounded function. 
The direct transform gives 


(p + x)B(p) = + ap + kj (f 2 + ff ' ) 


(33) 


and 


F(p) = 


g P 


(p+x) (p+y) 




a P + k 


p + x P + x 


m2 


(f + ff 1 ) 


(34) 


Applying the inverse transform to Equation (34) 


f(t) = + a e" xt 


y - x 


U- 

CL Lp + x 


G 


(f + ff) 


1 


f l + f 2 + f 3 + + + + 


(35) 
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where 

f (e xt - e •^ rt ) , -xt , 0 

f., = g — + ae . (3 

1 & y - x 

Equation (35) is more easily carried out by changing to Heaviside 
notation using the definition * 


j -' - , 

p + x 

and expanding the non-linear terms 


f l + f 2 + f 3 ^ + + = 


f x + kD 


[ f l 2 + f l f l’_ 


^kD|f 2 +2f 1 f 2 + f 1 'f 2 + f 1 f 2 * + f 2 fj 


hkD f 3 +2f 1 f 3 + 2f 2 f 8 


+ ' ^3 + f 3 f 2 ' + ^1^3 ' + ^2 f 3* + ^3^3 * 


+ + + + . 


In Equation (38) each bracket term on the right can be paired with 


a term on the left 


f 2 5 kD 


9 

f + f f 1 

z i Vi 


f 3 skD [ f 2 2 + 2f l f 2 + f l' f 2 

i 

+ f f 1 + f f ' 
12 T1 


f 4 skD f 3 ‘ + 2 tl f 3 + 2f 2 f 3 


+ f l' f 3 + f 3 f 2' - Vs’ +f 2 f 3' +f 4 f 3’ 
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Since f^ is known fg may be obtained by a "Heaviside shift" which 
may easily be programmed on a computer. When f ^ is known may be 
found. This procedure may be continued to obtain as many terms as 
required. 

The computer program is based on the following operation: 


/ pVy X 


, xt,n> 
(e t ) 


= f(t) 


( 42 ) 


where x and y are real and n is a positive integer including zero. 
There are two cases: 


Case I 


y = 

f(t) = 


-x 

xt ,n-l-l 
e t 

n + 1 


( 43 ) 


Case II y t -x 

-yt 

<~y>} <n+1) 

xt , n 

•e t 

(x+y)| n ni 

xt , n- 1 

-e t 

(x+y)^ (n-1) I 

-e xt t° 

{- (x+y) ^ (n+1) (0) J ( 

While the machine program handles only real values of x and y 
the results may be converted into sinusoidal quantities quite readily. 


f(t) = (ni) 


(- 
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6. THE PENDULUM 

The equation of motion for an idealized pendulum is 

— ■ -J* = - k sin 0 (4J 

dt 2 

where k is the ratio of the gravitational constant to the length of the 
pendulum. 0(t) is the angle of displacement in radians. 

14 

The solution for the period of oscillation is well known in this 

case and the expansion method would generate secular terms since sin 0 

would be approximated by the first few terms in 0. 

The derivative method also generates secular terms if the data is 

14 

sorted on k before summation. Knowing the period one may approxi- 
mate the rotational velocity as 




1 + <ir > 2 e 2 + <t^-> 2 ? 4 + < 4r!nr >" e 6 + + + 


where (3 2 = and b = 0’ (0). 


If data in p is then summed on the basis of the odd harmonics of 
Ul for f(0) = 0, f'(0) = b 


F(p) = 

(P +i ) 


_b 

8 


3 , 3b 5 , 3b 7 

512k 4096k 2 


P 

, 2 2„ 2 Q 2. 

(p +w 1 )(p +9u x ) 


15b 7 513b 9 

2048k 262,144k 2 


+ + + 


P 

(p 2 +u 1 2 ){p 2 +9u> 1 2 )(p 2 +25u 1 ‘ 


1,048, 576k‘ 


16,777,216k 


+ +• + 


(47) 
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The transforms may easily be evaluated to obtain the approxi- 
mation 


0(t) = — sin co -j t - 
0 )^ 1 


b 3 + 3b 5 
— + 


3b 


512k 4096k 2 


*•] 


• 3 * 
sm cj -^t 

6w 3 


G 


p; 7 q 

3b , 15b , 513b 

32 2048k 262,144k 2 


81b 


11 


1,048, .5 76k 3 
81b 13 

16,777,216k 4 


4* + +■ 


sin Wjt 


150u. 


(48) 

The maximum value of b = f ' (0) for which is valid is 2'vr. 

If this value is used for the series forming the coefficients they appear 
to converge very rapidly and it is assumed that they do converge. 

If sin Ujt is replaced by its maximum value one may then examine 
the first three terms of the alternating series. Again they appear to be 
absolutely convergent and the assumption is made that the result is ab- 
solutely and uniformly convergent and is therefore a solution for the 
original differential equation for the case where f(0) = 0 and f 1 (0) = b 
in radians per second. 

Figures (2) through (5) show analytic results compared with re- 
sults obtained by numerical integration for four cases. 





6(t) IN RADIANS 





6(t) IN RADIANS 




0(t) IN RADIANS 



to 

w 
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7. SECOND ORDER EQUATION 


One of the most interesting equations under study was the second 
order equation 

f" + (x+y) f* + xyf = kf 2 . (49) 

Here x and y are real and positive while k is a small positive or 
negative real number. 

In addition 

i m x / ny (50) 

where m and n are any positive integers. If this inequality does not hold 
multiple roots will appear causing secular terms in the final result. 
Evidently x and y must be irrational numbers. 

In the following operational transforms the initial conditions will 
be abbreviated as 

f(0) = a, f'(0) = b. (51) 

The first three terms of the complete solution by the derivative 
operational method are 
F(p) = 

p(ap+b+ax+ay) 

(p+x) (p+y) 

+ kp 
+k 2 p 


(p+x) (p+y) (p+2x) (p+x+y) (p+2y) (p+3x) (p+2x+y) (p+x+2y) (p+3y) 


(p a 2 +2pab+4abx+2b 2 +3pa 2 x+2a 2 x 2 +3pa 2 y+4aby+6a 2 xy+2a 2 y 2 ) 

(p+x) (p+y) (p+2x) (p+x+y) (p+2y) 

r* " 

AO o o on no nn nn 99 

2p a +16p a x+lOp a b+102p a xy+62p a bx+62p a by+20p ab + 
206pa 3 x 2 y+120pa 2 bx 2 +270pa 2 bxy+84pab 2 x+84pab 2 y+20pb 3 + 

132a 3 x 3 y+72a 2 bx 3 +252a 2 bx 2 y+72ab 2 x 2 +144ab 2 xy+24b 3 x+ 

_ 3 3 3 2 3 2 3 3 2 3 2 3 4 3 2 

24b y+16p a y+46p a x +56pa x +46p a y +24a x +206pa xy + 

2, 2 3 3 3 2 2 „„„ 2, 2 , 2 2 „„„ 3 3 

120pa by +56pa y +204a x y +252a by x+72ab y +132a xy + 

9 Q Q A 

L72a by +24a y 


( 52 ) 
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If this equation is evaluated it represents three terms of a formal 
solution of Equation (49). Each term multiplied by a given power of k is 
made up of negative exponentials added together. If x and y are irrational 
numbers there will be no repeated roots. If it is assumed that the infinite 
series of such terms behave in the same way some finite constant M m 
may always be chosen to replace each k multiplier and we have 


kM >f(t). 
n 


If these M n follow a normal pattern they eventually decrease with 
increasing n so that all M n may be replaced by the maximum value M. 


k n > f(t) 


Evidently the result, under these assumptions, converges abso- 
lutely and uniformly for t > 0 when k is sufficiently small. 

Fortunately it is not necessary to rely on so many assumptions 
since Poincares theorem assures convergence for sufficiently small 
values of k in any case. 

Now Equation (52) is the solution of Equation (49) in operational 
form. The operational solution is sometimes as interesting as the final 
evaluation for f(t). In this case the initial and final value theorems 
applied to the individual terms show that the final value of f(t) is zero 
for all three terms and the initial value is zero except for the first term 
where f(0) = a. 

A computer program evaluated Equation (52) directly for given 
values of x, y, f(0), f'(0) and t. The result is shown as Figure 6. 



1 


Numerical 

xxx Analytic 
(3 operational terms only) 


f ;/ +(x+y)f Vxyf-kf 2 
I I 

x = y =yir 

f(o) = l f (o) =0.5 

k«-o.i 


figure: 


8 10 IZ 14 16 18 Z0 22 24 26 28 

— seconds 
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Memory space in the computer is always at a premium with 
these solutions. While Equation (52) contained both initial conditions 
(a and b) the next term would have required more memory than available. 
One may obtain more terms (k n ) by sorting data either f(0) = a = 0 or 
f'(0) = b = 0. 

Equation (55) shows the - operational solution for f* (0) = b = 0 and 

f(0) = a 

F(p) = a(p 2 + xp + yp) + Ka 2 (p 3 + 3xp 2 + 2x 2 p + 3yp 2 + 6xyp ± 2y^p) 

(p + x)(p + y) (p + x)(p + y)(p + 2x)(p + 2y)(p + x + y)- 

2p-* + l6xp^ + U6x 2 p3 + I02xyp3 + 206x 2 ;yp 2 + 206xy^p^ + 132x 3 yp 

+ 20lpc 2 y 2 p + l6yp^ + U6y^p 3 + 56x 3 p 2 + 5>6y3p 2 + 2lpc^p 

+ H 2 a 3 b- 132xy 3 p + 2liy*p 

(p + x)(p + y) (p + 2x)(p + x + y)(p + 2y)(p + 3x)(p + 2x + y) 

(p + x + 2y)(p + 3y) 

10p 8 + X62xp^ + l62yp 3 + I090x 2 p^ + 394.2x 3 p^ + 13 ,U 66 x 2 yp^ 

+ 13,ii66xy 2 p^ + 39,922x3yp^ + 63 5 76l!X 2 y 2 p^ + lUbl36x3y 2 P 3 

+ lUU,136x 2 y 3 p 3 + 9l;lpc^y 2 p 2 + 219,?4b: 3 y 3 p 2 + 63,36Qx 3 y 2 p 

+ 118,65)6x^37^3 + 233 Cfeyp^ + I090y 2 p^ + 39U2y 3 p^ + 8260x^p^ 

+ 63 , 8 £ 6 x I *yp 3 + 39,922xy 3 p i; + 10,008x% 3 + £2,l£2x^yp 2 

+ 8260y^’p^ + 6U80x^p 2 + 16, 592x^yp + 63,856xy^ l p 3 + 10,008y^p 3 

+ 1728x' l? p + l31i,9U4x 2 y^p 2 + 52,15>2xy 3 p 2 + 6U80y^p 2 + 118 , 6 £ 6 x 3 y^p 

+ K 3 a^ . + 63,360x 2 y^p + l6,992xy^p + 1728y 3 p 

(p + x ) (p + y) (p‘ + 2x) (p + x + y) (p + 2y) (p + 3x) (p + 2x + y) 

(p + x + 2y)(p + 3y)(p + Ux)(p + 3x + y)(p + 2x + 2y)(p + x + 3y) 

(p + 4y) 
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8 Op 12 + 2'336xp i:1 - + 2336yp 11 + 30,288x 2 p 10 + 63,r?6xyp 10 
+ 7li9,8UQx 2 yp 9 + 7li9, SiiOxy^p 9 + 5,137,232x3yp 8 + 8,029,9Oljx 2 y 2 P 8 
+ 5,137 ,232xy 3 p 8 + 22,li67,296x^yp 7 + I 48 , 758 , l?6x 3 y 2 p 7 
+ i|.8,758,176x 2 y 3 p 7 + 22,li67,296xy^p 7 + 65,bi3,30ljx£yp^ 

+ l8i; J 7U9,l4UC!xV 2 P^ + 256,U59,56Gx 3 y 3 p 8 + l8ii,7lt9,iiUOx 2 y^p^ 

+ 128,382,OOQx 6 yp 3 + U52,36l,232xVp 3 + 8i6,660,896kMp 5 
+ 8l6,660,896x 3 y^p^ + U52,36l,232x 2 y 9 p^ + 7!li,697,l8ljx^y 2 pk 
+ 1,610, 209, 192x5^ + 2,093, 793, 712xW* + 1,610, 209,192x3^ 

+ 71ii,697 i l8bc 2 /p li + 701,61b,736x 7 y 2 p 3 + 1,917, 317, 088x 6 y 3 P 3 
+ 3,110,628,5l2x%V + 3 5 110 ,6285 5l2xV 3 P 3 + 1,917, 317, 088 xVp 3 
+ 388,007 ,OliOx 8 y 2 p 2 + 1,259, 06 l, 2 l 6 x 7 y 3 p 2 + 2^77,021,376x 6 A > 2 
+ 3 , 089,805, 696x 3 y 3 p 2 + 2,1*77,021,376 xV 6 P 2 + 92,0UU,800x 9 y 2 p 
+ 3 U 8 , 8 U 2 , 8 80 x 8 y 3 p + 8 lU , 798 , 080 x 7 y^p + 1 , 231 , 879 , 68 Qx 6 y^p 

+ 1, 231,879, 680x^y^p + 30,288y 2 p 10 + 229,92Qx 3 p 9 + l,13lj.,00Qx^p 8 
+ 229,920y3p 9 + 3,8lO,528x^p 7 + l,13U,000y^p 8 + 8,888,9liipc 6 p 6 
+ 65, 14i3,30lixy^p^ + 3,8l0,528y^p 7 + 8,888,9UUy^p^ + 167,2 38, 528x 7 yp^ 
+ lli,37U,2UQx 7 p 3 + 128,382,000xy^p^ + 138,276,288x 8 yp 3 
+ l5,769,920x 8 p^ + lU,37l4-,2U0y7p^ + 65, U73,92Qx 9 yp 2 + 11,166, 336x 9 p 3 
+ 167 ,238, 528xy 7 p^ + 15,769,920/p 1 * + 13,U78,liOOx 10 yp 
+ U,589,568x 10 p 2 + 7Ol,6I0i,736x 2 y 7 p 3 + 138,276, 288X/P 3 
+ ll,l66,336y 9 p 3 + 829,UUOx i:i p + I,259,06l,2l6x3y7 p 2 
+ 388 , 007 , OUOx 2 y 8 p 2 + 65,U73,920xyV + [ t ,589,568y 10 p 2 
+ 8ll;,798,080x^y 7 p + 3U8,8U2,88Gx 3 y 8 p + 92,0l4li,80Gx 2 y 9 p ^ 

+ l± 13,U78,l|Q0xy 10 p + 829,UU0y LI p 

(p + x) (p + y) (p + 2x) (p + x + y) (p + 2 y) (p + 3x) (p + 2x + y) (p + x + 2y) 

(p + 3y)(p + bc)(p + 3x + y)(p + 2x + 2y)(p + x + 3y)(p + liy)(p + 5x) 

(p + be + y)(p + 3x + 2y)(p + 2x + 3y)(p + x + l*y)(p + 5y) 


( 55 ) 
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Equation (56) is the operational solution for f(0) = a = 0 and 
f*( 0 ) = b. 

F(p) = bp + 2 Kb 2 p 

(p + x)(p + y) (p + x)(p + y)(p + 2 x)(p + x + y)(p + 2 y) 

+ K 2 b 3 ( 20 p 2 + 2 Upx + 2 lipy) 

(p + x)(p + y)(p + 2 x)(p + x + y) (p + 2 y)(p + 3 x)(p + 2 x + y)(p + x + 2 y)(p + 

600 p^ + 2760(xp3 + yp3) + 3888 (x 2 p 2 + y 2 p 2 ) + 8976xyp 2 

+ 1 + 662 l|(x 2 yp + x^p) + 1728 (x3p + y 3 p) ' 

(p + x)(p + y)(p + 2 x)(p + x + y)(p + 2 y)(p + 3 x)(p + 2 x + y) (p + x + 2 y) 

(p + 3y)(p + lpt)(p + 3 x + y)(p + 2 x + 2 y)(p + x + 3y)(p + Uyj 

'39 9 600p 7 + iiii3 s ?20(xp 6 + yp 6 ) + 1*972, 080 (xV + y*P^) 

+ U*27U*l6Gxyp? + U,l£7,376(x 3 pk + yV*) + 13*682* 128 (x 2 yp^ + xy 2 p^) 

+ 3*it03*l8lt(x^p 3 + y^p 3 ) + 27,li.02,336(x3yp 3 + xy 3 p3) 

+ 14)-*it28 ; , 70lpc 2 y 2 p3 + 3*3U3*U08(x^p 2 + y 3 p 2 ) 

+ 22,809*600(x^yp 2 + xjAp 2 ) + 33*2ll*8UO(x3y 2 p 2 + x 2 y3p 2 ) 

+ 7*237*600(x 3 yp + xy^p) + 22*732*000 (xVp + x 2 A>) 

+ K^b 3 _+ 32*6U7*680(x3y3p) + 829,lUiO(x^p 4- y°p) 

(p + x)(p + y)(p + 2x)(p f X + y)(p + 2y)(p + 3x)(p + 2x + y)(p + x + 2y) 

(p + 3y)(p + ipc)(p + 3x + y)(p + 2x + 2y)(p + x + 3y)(p + Uy)(p + 3x) 

(p + Iff + y) (p + 3x + 2y) (p + 2x + 3y) (p + x + Uy) (p + 5y) 

( 56 ) 

Figure (7) shows the results for the analytic solution compared with the 
results of numerical integration for k negative. Evidently the magnitude 
of k = -0. 3 is n too large” for convergence or more than 5 operational 
terms' are required. 




t seconds 


03 

© 
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The actual evaluation of Equation (49) as a time function is more 
quickly and more easily carried out by the expansion method, but in this 
case any numerical check solution must be programmed on the computer. 

As an example take the case where f(0) = 0 and f ‘ (0) = b, which is 
the evaluation of Equation (56). Notice particularly that the constant "a" 
is not f(0) as usual but a newly defined constant. 

Equation (49) may be written in Heaviside notation with one boundary 
condition as 

|p 2 + (x+y) p + xyj .f = bp + kf 2 (57) 

where f(0) = 0 and f * (0) = b 


f(t) = 


bp 

(p+x) (p+y) 


+ 


kf 

(p+x) (p+y) 


(58) 


1 

(p+x) (p+y) 

where 

1 

a = 

y-x 


(p+x) 


c 

p+y 


and c = 


1 

x-y 


f(t) =. abe *+ bee ^ + 



f 5 f 1 + f2 + f 3 + f 4 + + + 



(59) 


(60) 

(61) 


and 



+ (f 2 2 + 2f x f 2 ) 

+ < f 3 2 + 2t l f 3 + 2f 2 f 3> 

+ (f 4 2 + 2fjf 4 + 2f 2 f 4 + 2f 3 ( 4 ) 


H" + + 


(62) 
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f(t) " + f 2 + f 3 + f 4 + + and 


D a ["-X- + — r — ”1 (63) 

|_P+x p+y J 

f ^ = ab e + be e ^ (64) 

f 2 = Dk (f x 2 ) (65) 

f 3 = Dk (f 2 2 + 2f x f 2 ) (66) 

f 4 = Dk (f 3 2 + 2f 1 f 3 + 2f 2 f 3 ). (67) 


In carrying these operations out with the computer the two terms 

of are set up as shown in Equation (64). A multiply program then forms 
2 1 

f^- . All terms are then given a common denominator so they may be 

packed into the smallest number of terms. The next step is to multiply 

all terms by k. Each of the operators of Equation (63) is then applied 

2 * 

individually to kf ^ . A common denominator is then found- so that the 
terms may be packed. The result is then read out in decimal and punched 
as a hexadecimal tape for future operations. 

The results are carried in the machine as algebraic quantities and 
exponentials so there is no need for laborious hand calculations. 

In this particular case no secular terms appear since it is assumed 
that the inequality of Equation (50) is satisfied. 
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3 

The approximation carried out through the k term only is 
given by 

f(t) = f 1 + f 2 + f 3 + f 4 + + + . 

Where 

= ab e"' x *' + be e*"^ 



- aV^e-^-l) - 2a 2 C e“ xt (e-y t -l) + acV** (e (x_2y)t -l) 

x y x-2y 

- ^ce-ytCe^^-l) - 2ac 2 e-y t (e- xt -l) - cV^ (e^-l) 

, X y 



k 2 b 3 


+ aV^e^-l) - 2a 5 e- x ^(e~ xt -l) + ha^ce"^ (e~ (x4y)t -l) 

x 2 x 2 y(x+y) 

- Ua^ce -xfc ( e -xt -1 ) - aVe^e** 2 ^-!) + 2a3 c 2 e" xt (e’- xt -l) 

xy y(x-2y ) x(x-2y) 

+ a^ce~’ x ^(e” 2 x ^-l) - 2 a^ce"‘ x 1 : ’(e” 5 r t'-l) + I*a 3 c 2 e” x ^(e* , ‘( x+y ^-l) 

x ( 2x-y ) y(2x-y) x(x+y) 

- ita3c 2 e~ xt ( e -1 ) + a 2 c^e“ x *'Ce“ 2 ^ r * : '-l) - 2a 2 c3e -x ^(e”y^-l) 

xy y 2 y 2 ' 


+ 2 aW xt (e-( x+ y) t -l) - 2 a li ce" xt (e~y t -l) + 2 a 3 c 2 e _xfc (e- 2 y t -l) 

xlx+y) xy y2 

- Ua 3 c 2 e"* t (e~y fc -:i) + 2a 2 c 3 e~ xt (e( x ~ 3y - >t -l) + 2a 2 c 3 e~ xt (e^-l) 

" (x-2y)(x~3y) y(x-2y) 

+ 2 a 3 c 2 e~ x ^(e~( x+ y)^-l) + 2 a 3 c 2 e" x ^(e^ x “ 2 y ^-l) 
i2x-y J (x+y) (2x-yJ (x-2y) 

+ 2a 2 c 3 e ^(e 2y ^-l) + iia 2 c 3 e _x ^(e^ x-2y ^-l) 
xy xlx-2y) 

- 2 ac^e" xfc (e^ x “ 3 y) t -l) + 2 ac^e- xfc (e (x " 2 y 5 t -l) 

y(x~3y) y(x-2y) 


Cont. 
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Continued 


+ kb 


-2v3 


+ kV 1 


+ 2a^ce-y fc (e^~3x)t^ 1) _ 2 a ) ^c.e-y t 

x(3x-y) x(2x-y) 

+ 2a 3 c 2 e“ yt (e“ 2xt -l) - ItaVe^fe^H.i) 
xy y(2x~y) 

- 2a 2 c 3 e" yfc (e-( x ‘ f y)t_ 1 ) + 2a 2 c 3 e”y t (e(y- 2x )t-l) 

(x-2y ) (x+y ) (2x-y ) (x-2y) 

+ 2a 3 c 2 e“ yt (e(y- 3x ) t -l) - 2a 3 c 2 e”y t (e- xt -l) + 2a 2 c3 e ”^(e“ 2xt -l) 
(3x-y)(2x-y) x"(2x-y) x 2 

- Iia 2 c3e“3 r * : ’(e“ x ^-l) + 2ac^e~^ (e"^ x+ y^-I) - 2ac^e~‘^ (e** 3 ^-!) 

x 2 y(x+y) xy 

+ a^e-^Ce- 2 ^-.!) - 2a 3 c 2 e~y t (e^-l) + ^cV^fe-^ 4 ^-!) 
x 2 x 2 y(x+y) 

- Ua^c 3 e~y^(e -x ^-l) - ac^e^^Ce -2 ^-!) + 2ac^e’“y^ , (e” x ^-l) 

xy "Tl^Tx^y} x(x-2y) 

+ a 2 c 3 e~yk(e” 2:jct -l) - 2a 2 c 3 e~y t (e”^-!) + Uac^e” y ^(e“'^ x+y ^-l) 
x(2x-y) ' y(2x-y) ' x(x+y) 

- Uac^-e ( e*~y^ -1 ) + - 2c^e“y^(e~y fc -l) ‘ 

- xy y 2 y 2 *^ 

- a V^Ce- 3 **-!) + aV^e^-i) - Ua^e^e^ 2 ^*-!) 

3x 3 x 3 xy(2x+y) 

+ 2a 6 ce“ xfc (e“ 2 * t -l) + 2aVe“ xt (e-fr^H-P - a^c 2 e~ xt (e~ 2xt -l) 
x^y x(x**2y)(x+2y) x 2 (x-2y) 

- 2a 6 ce- xt (e- 3xt ~l) + 2a^e~ xt (e^^H-l) 

3x 2 (2x-y) x(2x-y)(x+y) 

- ^cV^Ce^ 2 ^^-!) + Sa^e^Ce"^^-!) 

x 2 (2x+y) ' x 2 (x+y) 

- 10a ii c 3 e“ xt (e-( x+2 y)t-l) + lOa^cV^ (e-k^-l) 

xy(x+2y) xy(x+y) 


Cont . 
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Continued 
+ k 3 ^ 


- a 7 e~^(e~ xt -l) + - l ia 6 c e " xt (a'* t -l) 

x3 . xy(x+y) x^y 

- ^e^ je^-l) + ZaVe^fe^-l) + a^e^ (e“ 2xt -l) 

xy(x-2y) ” x 2 (x-2y) x 2 (2x-y) 

~ 2a 6 ce~ xt (e~y t -l) - ka^cV^e"^-!) + (e- 2 ^-!) 

xy(2x-y) x^y xy 2 

- lOa^cV ^Ce-yfe-p - Ua^c 2 e' xt (e“ (x+2 y) t -l) 

xy 2 y^(x+2y) 

+ 8a^c 2 e~* x ^(e**^ x+ y^~l) + Ua 1 ^ V xt (e^-l) 

y^Oc+y) Jy 2 (x-2y) ~~~ 

- uA3 e -xt( e -(x4y)t < _i ) .i )a 5 c 2 e -xt (e -(2xtr)t„i ) 

y(x-£yj$c+y) “ y(2x-y) (2x+y) 

+ iiaVe^Ce- 2 ^-! ) - (e“3yt_i) + l ia 3 c li e“ xfc '(e’ 2 y t -l) 

2y 2 (2x-y) ^3 

- UaVe^e^-l) - 2a^c 3 _e“* xt (e“ 2 y t -l) + uA^e" 3 * (e _xfc -l) 

xy 2 y 2 (x-2y) xy(x-2y) 

+ 2a^o 2 e” xt (e" 2xt ~l) - uVe'^ (e -yt -l) - Ua^e"^ (e-y^l) 
xy(2x-y) y 2 (2x-y) ^3 

+ a3 c ^e~ xt (e^y> t -l) + a^V^e" 2 ^-!) 

(x-Uy) (x-2y)2 y(x-2y) 2 

+ 2a^c 3 e~ xt ( e -(x+2y)t^ l) + 2a^c 3 e _xt (e^^^-l 
(x-2y) (x+2y ) { 2x-y ) J2x-y)(x- 2y){x-3y) 

+ U 3 cU e -xt( e - 3 yt^i) + ^c^e^Ce^” 3 ^-!) 

3xy(x-2yJ x(x-2y)(x-3y) 

- 2a^cV' xt (e^ x " J+ ^ t -l) + 2a 2 c 5 e~ xfc (e( x “3y)t_. 1 ) 

y(x-UyHx-2y) y(x-3yMx-2y) 4 

- a 3 c^e~ xtl (e~ xtl -l) - a^c^e"^ (e“ 2xt -l) + 2a^c 3 e"’ xfc (e~ yt -l) 

x(x-2y) 2 x(2x-y)(x-2y) y(2x-y) (x-2y) 


Cont. 
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f 


3 


Continued 


+ k 3 b^ 


- U 3 c^e~ xt (e^ x4 y) t ~l) + Ua^e^e-^Ce^-l) - aVe^Ce’ 2 ^-!) 

x(x-2y)(x+y) xy(x-2y) y 2 (x-2y) 

+ 2a 2 c 5 e- xt (e-y t -l) - + 2 a 5 C 2 e- xt (e"( x '*y> t -l) ' 

y 2 (x-2y) " 3x(2x-y) 2 (2x-y) 2 (x+y) 

x(2x-y} (2x-fy) x(2x-yHx+y) 

- 2a3A“ xfc (e“^ +2 y) t -l) + 2a3 c ^e- xfc (e“( x+ y^-l) 

y(‘2x-y) {x+2'y) " y(2x-y)7By) 

+ a^cV^efr- 2 ^-!) + 2a^c 3 e~ ac ^(e‘“^^~l) 

(2x-y) 2 (x-2y) xy(2x-y) 

+ lia^c3 e - xt (e( x - 2 y) t -l) - 2a3A- xt ( e ( x ’ 3 y) t >l) 
x(2x-y ) fx~2y; yC2x-y)'(x-3y) 

+ 2a3 c ^- xt (e^- 2 y )t -l) - i; a 3 c U -xt (e -(x+2y)t , 1} 

y(2x-y) (x-2y) x 2 (x+2yY~ 

+ U 3 c i4 e- xt ( e - 2 y t -l) ~ W^eJ^e" 3 ^-!) + 2a 2 cV^e- 2 ^-!) 
x'T 3xy 2 ~ xy 2 

+ Ua 3 c^e~ xfc (e Cx - 2 y^ t -l) - Ua 2 c g e“ act (e^- 3 y)t_ 1 ) 

x 2 (x-2y) xy(x-3y) 

+ ^o^e^Ce^yH-l) + ac 6 e- xt (e^y) t -l) 

xy(x-2y) y 2 (x~Uy) 

- 2ac 6 e- x fc(e^ x -3y)t_ 1 ) + &c 6 e -xt (e^ 2 ^ -l) ■ 

y 2 (x-3y) y^(x-2y) 

- a^e-^Ce^y-^^-I) + 2a 6 ce-y t (e(y3x)t,i) 

x 2 (Ux-y) x 2 (3x-y) ' 

- i^a g c 2 e~y t (e~3xt ,p + i^cVy^ety^t _p 

~ 3x 2 y xy(3x-y) 

+ 2 a^c 3 e“y t (e“^ 2 x+ y) t ~l) - 2 a^c 3 e“y b (e^ y " 3 x ^ t -l) 

x(2x+y)(x-2y) x(3x-y)(x-2y) 

- 2 a^c 2 e-y t (e(y ijx ^~l) + a ^e^( e ^ 4) - j^e-yt^xt,!) 

x(lpt-y)(2x-y) x 2 (2x-y) 3x3 


Cont. 
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3 


Continued 


+ 


+ W+cVyt(e- 2xfc -l) - 10a 3 c^e“-y t (e~^ 2x ' f y> t -l) + 5a3 e (e^-l 
x3 xy(2x-y) x 2 y 

~ a 6 ce-yt(e(y” 2x H.i) + 2 a *c 2 e^ fc (e- gxfc -l) 
x 2 (2x-y) x 2 y 

- Ua 5 c 2 e-y t (e ( y~ 2x)t -l) - 2a 1 *c3 e -y t 

xy(2x-y) x(x-2y)(x+y) 

+ 2a k C Vyt(efr- 2x > t -l> + 2a ? c 2 e“y t (e^-3x)t_ 1 ) 
x(2x-y ) (x~2y) — xT3x-y) (2x-y)- 

- 2a^c 2 e-y t (e- xt ~l) - Ua^e"^ (e^-l) + 

x 2 (2x-y) xy(x-^y) 

- lOa^c^e**^ (e -3 ^-!) - Ua^c^e“y^(e“^ 2x+ ^ r ^-l) + Ua^c^e"’^' (e _2x ^~l) 

x^y y 2 (2x+y) xy 2 

+ ^ 3 c lf e~y t (e~^ +2 y)t-l) - 2 a 3 c li e“ yt (e~ ,2xt -l) - kahch^Ce^-l) 

y(x-2y) 2 xy(x-2y) 3xy(2x-y) 

y(2x-y)(x+y) y 2 (x+2y) 

+ Sa^^e**^ (e~( x ^y )^-l) ~ Ua^c3e“yt(e^3 r "“ 2x )^'— 1) 

^(x+y) y 2 (2x-y) 

- Ua 3 c^e’ , y fc (e“( 3ct y) t -1) + Ua^A-^Ce^-^^-l) 

y(x-2y)(x+y) y(2x-y ) (x~2y) 

+ UaMe-^Cefr^xH-i) . - haVe"^ (e** 0 *-!) 

y(3x-y)(2x-y) xy(2x-y) ^2 

- a 2 c g e~yt(e-3yt,l) + 2a 2 c^e"yt( e -( x ty)t_ 1 ) 

3y(x~2y) 2 (x-2y) 2 (x+y) 

+ 2a 3 c^e-y t (e-( 2x4 y ) t ~l) - aSc^e-^Ce" 2 ^-!) 

(2x+y)(2x-y)(x-2y) y(2x-y) (x-2y) 

+ Ua 2 c^e-’^(e“( x+2 y) t -I) - 2a 2 C 3 e -yt(e~ 2 yt.i) + 2ac 6 e“^(e-3yt>i) 
x(x-2y)(x+2y) xy(x-2y) 3y* (x-2y) 

- ac^e*”y^ (e“ 2 y^-i) - a 2 c^e _ y^(e^y" 2x ^-l) 

y 2 (x-2y) (2x-y) (x=2y) 2 


Cont. 
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f 3 


Continued 


+ kV* 


- 2a3 c ^e-yt( e (y-3x)t_ 1 ) + 2a3 c U e -y t (e~ xt ~l) - 2a 2 c g e~y t (e~ 2xfc ~l) 

(3x-y)(2x-y)(x-2y) x(2x-y)(x-2y) x^(x-2y) 

+ l ;a 2 c g e~y t (e~^-.l) - 2ac 6 e“y t (e"^ x+ yH, 1 ) + 2ac 6 e“y t (e’* cfc -l) 
x 2 (x-2y) y(x-2y)(x+y) xy(x-2y) 

- a^ 3 e -yt( e (y-bx)t_p + a^cV^e- 2 **-!) - 

(ipc-y) (2x-y) 2 x(2x-y)^ 3x 2 (2x-y) ~~~ 

+ 2 a 3 A~y t (e“ 23 ct -l) - 2 a 2 c^e~y t (e-( 23 d y) t -l) + a^V^Ce" 2 *^!) 
x 2 (2x-y ) y(2x-y)(2x+y) xy(x-2y) 

- a It c 3 e"y t (e~y^-i) + - iia 3 ^"^ (e^-l) 

y(2x-y) 2 x(2x-y)(x+y) xy(2x-y) 

+ a 2 c^e-y t (e~ 2 y t -l) - 2a 2 cVy t (e’-y t -l) - Ua 2 t£>-y* ( e -( 2 *+yH-l) 

y^(2x-y) y2(2x-y) x 2 (2xty) 

+ 8 a 2 c g e"y t (e-^ + y) t -.l^ - UaA^Ce’^y^-.!) 

x 2 (x+y) xy(x+2y) 

+ Uac 6 e‘*y^ (e**^ x+ y )^~1) - I^c^e-yt (e^-l) + uA^e" 2 ^-!) 
xy(x+y) ' x^y 2xy^ 

- iiac 6 e-yb(e“y^-l) - cV^ (e^yt.p + c 7 e -yt (e -2yfc._p 




f 4 is shown in the Appendix. 

For given numerical values of x and y one could probably show 

uniform convergence for this result. It is much easier to depend upon 
// 

Poincares theorem for k small enough to assure convergence. 
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8. ASSIMILATION OF SECULAR TERMS 


Transform solutions are formal solutions of the differential 
equation. When repeated roots occur in the operational transform 
secular terms appear in the solution. In cases similar to the pendulum 
where the period (for oscillatory functions) is known the summation of 
p n may be carried out by the computer, or if a problem is first solved 
by hand the summation can be carried out by the machine but there would 
be no point in wasting machine time. 

One may adopt the point of view that with machine solutions a 
sufficient number of secular terms may easily be obtained to approximate 
the solution or one may attempt to recognize the secular terms in such 
a way as to assimilate them into the solution as known functions. As an 
example of the two methods consider the equation 

f" + u 2 f = -kf 2 . (69) 


A formal solution to this equation may be obtained by the derivative 
process or by the expansion method. 

The derivative process, after summation and sorting gives Equation 
(70) directly while the expansion method (after appropriate changes in 
constants) gives Equation (71) directly. 


F(p) = 


ap 


2 ,. 2 

p +co 


ka 2 (p 2 +2(p 2 ) 
(p 2 +w 2 ) (p 2 +4w 2 ) 


9^4 99 A 

, k a (2p +10u p - 12u ) 

(p +io ) (p +4w ) (p +9u ) 

k 3 a 4 (10p 6 +150(o 2 p 4 +440id 4 p 2 +1200to 6 

, 2 2,2 T 2 , “ 2,2 , 2 “ 2 , . 2 , ie 2 , 

(p +co ) (p +4to ) (p +9w ) (p +16(o ) 

+k 4 a 5 (80p 12 +3,880io 2 p 10 +65,040(o 4 p 8 +444,440io 6 p 6 

+942,880(o 8 p 4 -3,040,320(o 10 p 2 -5,760,000io 12 ) 

, 2 2,3 . 2 2,2 . 2 Q 2,2 . 2^ 1C 2,2 . 2^ 2, 

(p +(0 ) (p +4u ) (p +9to ) (p +16io ) (p +25u ) 


+ + + 


( 70 ) 
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When F{p) is evaluated as a time function: 


f(t) = a cos iot 
2 


+ 


ka 

6w 5 


(cos 2iot + 2 cos cot - 3) 


,2 3 

+ - a 4 -- (3 cos 3 tot + 16 cos 2 ut + 29 cos tot + 60 tot (sin tot) - 48) 
144io 


, 3 4 
k a 


(5 cos 4 tot + 45 cos 3iot + 480 cos 2iot + 595 cos tot 


2160io + 3 o0iot (sin 2 tot) + 900 tot (sin tot) - 1125) 


(71) 


This equation may he rewritten as 
f(t) = 

2 


(a + 


+ ( 


ka 

3to" 

ka' 

6to' 


+ 2 ; i: : y + + + > cos u t 


+ 


144co 

2 
c 

144u? 


1fil 2 3 
16k a 


+ ( jjay + 


2160to 6 

+ 48 - ^ g- + + ) cos 2 tot 
2160io 

-t- + + ) cos 3iot 


144io‘ 


2l60io 


c t 3 4 

+ ( + + + ) cos 4tot 

2160(0° 

T v 2 3 lA 4 i 

+ (60wt) + (900iot) - 


L 


144to" 


2160to 


sin tot 


(300iot) k 3 a* . Q . 
+ _i ^ . S in 2iot 

2160 io 


ka 


2io 


48k 2 a 3 , 1125k 3 a 4 ,, , 

+ 4 — + ~ ~ : ~s — 1 


144w 


2160u) 


(72) 
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To assimilate the secular terms (sin wt and sin 2wt) the relation 


cos(w-i-p) t = cos tot cos pfc - sin wt sin j3t 

. 2.2 A A 


= cos wt (1 


_P_L 

2 ' 


(3 t 


= sin wt ( 


Pt 




4 1 
+ + ) 


+ + + ) 


1 1 3 ' 

is compared with the cos wt and sin wt terms of Equation (72). 


(73) 


(a + 


+ t 


ka 


3w 


<60w) 


^ 29k 2 a 3 , , , 

+ -r — + + ) cos wt 

144w 


k 2 a 3 (900w) k 3 a 4 

a + ~a — 

144w 2160w 


sin wt 


(74) 


Divide the coefficient of (t sin wt) by the coefficient of cos wt. The 
result is; 


K , 2 2 _.33 

5 k a . 5 k a 

IS *• rr ~ 


12w‘ 


18w" 


(75) 


therefore in Equation (73) p is the negative of this value. 

The new rotational velocity is 
K1 2 2 .,3 3 

, a ,, 5k a 5k a , , . /r __. 

W + P = w (1 - 2 — “ e — + + ) - w. • (76) 

12w 18w 

The coefficient of (t sin 2wt) when divided by the leading coefficient 
of cos 2wt gives 

2 

\ 12w 4 / 

which matches the second term on the right of (76) since this is double 
frequency. 
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The period is 


T = 



(77) 


Figure 8 shows the graphical result of the solution by numerical integration 
for k = 0. 1, co = 1 and f(0) = 1. 


T 

( 1) (1 - 


2ir 

5(0. 01) 
12 


5(0.001) 

18 


6.283 
1 - 0. 0045 


= 6.31 seconds 


which appears to check the period in the figure. 

In Figure 9 f(0) has been changed to 3. In this case 


T*= o— =6,58 seconds 

5(0.01) (3) 5(0.001) (3) 

12 ~ - 18 

which is probably very close to the period of the numerical solution 
for this case. 

The solid line therefore represents the solution with the secular 
terms assimilated for both figures. 

Equation (71) may therefore be rewritten dropping the sine terms 
and changing to to to ^ in the cosine terms. The to terms in the denominators 
are unchanged. 

ka 2 

f(t) = 2 ~ (cos 2io^ t + 2 cos Ujt - 3) 

6u 

k 2 3 

+ ... (3 cos 3u).,t + 16 cos 2u,t + 29 cos u,t - 48) 

144w ill 




7 3 4 

k a 


2160io 


(5 cos 4(o ^t + 45 cos 3u^t + 480 cos 2io^t 
+ 595 cos io^t - 1125) 


(78) 


The crosses in both Figure 8 and Figure 9 indicate the equivalent 
solutions for Equation (71) including the secular terms as shown in this 
equation. 

2 

Convergence is assured if the factor ka/u is small enough. 
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9. THE VAN DER POL EQUATION 


The van der Pol equation is perhaps the most interesting equation 
studied. 

f" - A( 1 - f 2 ) f' + Bf = 0 (79) 

The work on the computer was actually carried out with two 
different equations 

f" - Af 1 + Bf = -kf 2 f * (80) 

where k = A in the solution 
and 

f" - 2Xf • + (X 2 - Y 2 ) f = 0 81) 


where 


X = 


A - f 2 k 
2 


Y = 


(A - f k) 


- B 


82) 


X and Y are time functions. 
In the final results 


A. “ a k , 

X(0) = = x and 


Y (0) = J (A ‘ a- k > - B = y 


where a = f(0) and b = f f (0). Again k = A for the van der Pol Equation 
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The first four operational terms for Equation (80) for the case 
f(0) = 0 and f T (0) = b are listed as Equation (83). 

F(p) = bp - 2Kb 3 p 2 

(p 2 - Ap + B) (p 2 - Ap + B)(p 2 - 3Ap + B + 2A 2 )(p 2 - 3Ap + B) 

+ 60K 2 b^p5 + 3U8K 2 b^Bp 3 - 2g2K 2 b^Ap k - U80K 2 b^ABp 2 + 2U0K 2 b^A 2 p 3 

(p 2 - Ap + B)(p 2 - 3Ap + B + 2A 2 )(p 2 - 3Ap + 9B)(p 2 - 5Ap + B + 6A 2 ) 

(p 2 - 5Ap + 25B) (p 2 - 5Ap + 9B + UA 2 ) 


bSliOK^b^p 10 - 112,680K 3 b7Ap9 + 7h3 3 352K 3 b 7 A 2 p 8 + 209,l60K 3 b 7 Bp 8 
2,5l2,536K 3 b 7 A 3 p 7 - 2 ,U93 5 614.8 K 8 b 7 ABp 7 + ui587,02UK3b 7 A> 6 
+ 11, 329j032K 3 b 7 A 2 Bp 8 + l 3 6l4-6 5 280K 3 b 7 B 2 p^ - l;,290 3 2UOK 3 b 7 A^p5 

- 2U s U11,26UkVa3b p ^ - 13,073,832 kVab 2 p^ + 1,612 ,800K 3 b7A 6 pk 

+ 2U,876,ii80K3b7A^Bp^ + 35,5l7,072K 3 b 7 A 2 B 2 p ii + U,005,720K3b7B3pk 
9,676, 800K3b7A^Bp 3 - 38,U25,920K 3 b 7 A3B 2 p 3 - lit , 96U , U 80K 3 b ? AB 3p3 
1+ Ul,112,OOOK 3 b7A i| B g p 2 '+ 12,868, 800K 3 b7A 2 B 3 p 2 - 756,000K 3 b 7 B 1 +p 2 ' ' j 


(p 2 - Ap + B)(p 2 - 3Ap + B + 2A 2 )(p 2 - 3Ap + 9B)(p 2 - 5Ap + B + 6A 2 ) 
(p 2 - 5Ap + 25B)(p 2 - 5Ap + 9B + UA 2 ) (p 2 - 7Ap + B + 12A 2 ) 

(p 2 - 7Ap + 25B + 6A 2 )(p 2 - 7Ap + U9B)(p 2 - 7Ap + 9B + 10A 2 ) 


(83) 


This was the first second order equation attempted and no 
s umm ation was obtained for the general case where both initial condi 
tions were present or for the case f(0) = a and f 1 (0) = 0. 
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Figure 10 shows the analytic solution for the case f(0) = 0, 
f'(0) = 1, K = 2. 25, B = 1 and A = -2.25. (Since A is negative the equation 
is not in the standard form for the van der Pol equation, ) The solid 
lines, as usual, are the numerical solutions. The agreement is surprisingly 
good since both numerical and analytic results are approximations. 



* Numerical solution 

x x x Analytic solution 
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The first three terms of the operational solution for Equation 
(81) for f(0) = a and f'(0) = 0 are 

F(a,p) - - 2 *P> 

Mn 

P* 6x^> 3 + i2 x 2 y 2 p3 _ 6y^p3.+ _ 88x 3 y 2 p 2 

Lt- l*lpty% 2 - $kx 6 p + llipcV 2 p - 66x 2 y\j + 6y 6 p_ 

(Mu) 2 M 3 iM33 


where 


+K 2 a^ 


~?52x 6 p 12 + 7?6xV 2 p 12 - 7£6xV* p 1? + 252y 6 p 12 + 10,360x 7 p 11 
-31,080x^y 2 p i;i + 31,080x3y^p 1 l - 10, 360xy 6 p 11 - 183,11*8x^3-0 
+ 553,O08x 8 y 2 p3-0 - 560, 136xVV° + 193 ,8l*0x 2 y 6 p 10 - 3561*y 8 p3-° 

+ l,835,Ol*Ox 9 p 9 - 5,623,l*88x 7 y 2 p 9 + 5,860,22l^xV l P 5, - 2 } 190,li { i l x 3 y 6 p9 
+ ll8,368xy 8 p 9 - ll,588,l*72x 10 p 8 + 36,1*30, 81*0x 8 y 2 p 8 - 39,775,728x 6 y it P 8 
+ 16,626, 86 JjxV 6 P 8 -l,707,5lilpc 2 y 8 p 8 +i* 8 ,Ul 6 ,I* 6 l*x :U T ? 7 -l58,060,0l*8x 9 y 2 p 7 
+ l8l*,083,360xV*p 7 -88,031i,l;32x5y 6 p 7 +H*,Ol 6 , 656 x 3 y 8 p 7 -1*02 , 000xy 10 p 7 

- 136,102,200x 12 p 6 + li67,291,696x 10 y 2 p 6 - 589,682,02l*x 8 yV’ 

+ 326,3hl,l52x 6 y 6 p 6 -72 3 3l3 } 272xV 8 p 6 + l i ,i l 89,90i|X 2 y 10 p 6 + 253 3 l6O,332x 1 V 

- 932,6l?,728x 11 yV + 1,290,537, 9 8Ux 9 yV - 828,06l,l*l*0x 7 y 6 pS 

+ 239,608, iil6x^y 8 p ? -25,137,920x 3 y 1 V +i*90,336xy 1 V -306,i*09,6l2x :i V* 

+ l,20l 1 ,893,7i|8x 12 y 2 p^ -l,8U2,131,132x 10 yV* + !,36l*,l*10,308x8y6pU 
-ii93 , 25U , 8 20x 6 y 8 p^+ 75,809, 561jx V 10 P^~3 ,335, l*l*i*x 2 y 12 pk+ 213 , 1*91 , 83 2x3£ p 3 

- 910,3ii9,256x 13 y 2 p 3 4 - l,5Ul,77U,0lj0x 13 y^p 3 - 1, 307,091*, 50l*x 9 y 6 p 3 
+ 571,950,696x 7 y 8 p3-ii9,i]5 } o96x% 10 p3 + 9,538,632x3yl2 p 3„65,678,9l*Ox l 6 p2 
+ 3 Oli, 91 7 , 120x 1 ^ l y 2 p 2 - 570,702,960x 12 y^p 2 + 51*7,560,000x 10 y 6 p 2 

- 282,79li,760x 8 y 8 p 2 )-7?,559,680x 6 y 10 p2 -9,?l5,280xV 12 p 2 +ll*,0l*0yl0p8 

- 23,256y 12 P 6 +17,388y 1 ^p^ -196,3]*axy 1 Jl p 3 + 360,000xW -l* 860 y l6 P 2 

y 2~~ '~ ' 2 — 

( M ll) (M 31 ) (M 33 ) 


M wn — (p 2 - 2wxp + w 2 x 2 - n 2 y 2 ) 


(84) 
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and the first four terms for f(0) = a and f'(0) = b are 


F(b,p) = bp -■ 2Kb 3 p 2 

M 11 M ii M 3i M 33 


+ K^b^ (60p^ - gOlgpk + I308x 2 p 3 ~ 3U8y 2 p 3 - 960x 3 p 2 + 96Qxy 2 p 2 ) 

%l M 31 M 33 M 5l M 53 M 55 


+ kV 


681;Op 10 + 225,36Qxp? - 3,l82,563x 2 p 8 + 209,l60y2p 8 + 25,o87,58lpt 3 | 
- h, 987,296xy 2 p^ - 120,35U,792xkp 8 + U8,6o8,688x 2 y2p^ - l,6U6,280y^]: 
+ 358,725,U56xV* - 2k7,585,liUOx 3 y 2 p 3 + 26 ,lU7 ,66i|xyV^ 

~ 6U7,316,888x 8 p^ + 69lt,177,Ul6xky 2 p^ - l5U,085,UU8x 2 y^)^ 

+ 1^,005,720^?^ + 6U6,993j920x^p 3 - l,0lU,259,200x 3 y 2 p 3 
+ 397>19i;,2iiQx 3 y^p 3 - 29,928, 960x/p 3 - 2 76,511, 200x 8 p 2 
+ 602,985,60Cbc 8 y 2 p 2 - 375,68l,600x^y^p 2 + U8,)j5l,20Gx 2 y6p“ 

+ 756,000y 8 p 2 ” ' 

M ll M 31 M 33 M 5l M 53 M 5^ M 71 M 73 M 75 M 77 


Here 

= (p 2 - 2v?xp + w^x 2 - n 2 y 2 ) (85) 


No satisfactory evaluation of Equation (84) has been obtained due 
to the large number of multiple roots. A computer program is available 
to evaluate this equation as a time function for given numerical values of 
the constants and the initial condition but unfortunately there appear to 
be too few secular terms to approximate the function. 
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The expansion type of solution might be very helpful here in 
identifying the series of secular terms and attempting to assimilate 
them. 

Equations (84) and (85) were generated with the idea of simplifying 
the root structure of Equation (83). Notice that the roots for Equation 
(84) are 

, 2 „ 2 2 2 2 . 

m um " 2toxp + u x - n y ) 

= |p 2 - « (A - a 2 k)p 

2 2 

+ (u 2 - n 2 ) ( A 2 a k ) + Bj . (86) 

where w and n take on the values shown in (84). Since a = f(0) this means 
that at least one initial value occurs in the root structure. 

Equation (85) is the same as Equation (83) if the values for x and 
y are substituted in the equation f(0) = a = 0. 

Poincare^s theorem insures convergence foi* small values of k 
so the result is an analytic function. 
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1 0 . DUFFING ’ S EQUATION 


The general form of Duffing’ s Equation is usually given as; 

f” + cf 1 + rf = -.kf 3 + G cos wt. (87) 

If the general solution for this equation was of sufficient impor- 
tance it could no doubt be obtained. In Heaviside notation with two initial 
conditions the original equation becomes 

(p+x) (p+y) f = ap 2 + bp + (x+y) ap + G cos ut - kf 3 (88) 


where a = f (0) and b = f’(0). 

2 . , / .... . v _ 


(p+x) (p+y) 


G cos wt 

1 

k / 

(p+x) (p+y) 

(p+x) (p+y) 



. (89) 


Even an approximate solution through the k term for this equation 
would take a great deal of time on a large computer. To reduce com- 
puter time take a = b = 0 and y = x so the equation now becomes; 


f = 


G cos wt 

2 . 2 
p + x 


k f 


2 2 
p + x 


(90) 


The first term on the right hand side is 
G (cos tot - cos xt) 


f i- 


V 


f 3 = 


, 2 2 \ 

(x - u ) 


- k 


- k 


— — 



1 


f 3 

2 2 


h 

p + X 




r 


2 2 
p + x 


3f l 2f 2 + Sf l f 2 2 + f 2 2 


(91) 


(92) 


(93) 


f - f^ + f2 + fg + + + • 


(94) 
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Computer results give: 



1 

4x 


cos 3xt - 


18x 



cos ut 


2x 

(x 2 - 9to 2 ) 


cos 3ut + 


3x 

2w (x-w) 


cos (x-2w)t 
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Rearranging the terms in (91), (95) and (96); 


f = 


G cos iot 9 kG 


kG. cos 3iot 


(x 2 - u 2 ) 4 (x^-io^) 


2 TT~ COS wt " ; ; 2 2,3 . 2 „ 2 


G 


, 2 2 . 
(x -to ) 


cos xt 


H 


4 (x _-io ) (x -9co ) 
2 


9 k G 2 t 

„ , 2 2.2 , 

8x (x -to ) 


) 4t 


+ + 


. . . 9k G 2 t , , 

sm xt ~ — ij— q— . 4- + 4- 

• 8x (x -to ) 


k G‘ 


n Q 2 , 2 2.3 

32x (x -to ) 


cos 3xt (1 + +) 


- 3^ 9 k ? K * + + + ) sin 3xt 


^8x (x 2 -u 2 ) 2 
3 k G 3 


sin 3xt 


16w (x-to) (x 2 -w 2 ) 3 


cos (x - 2w )t (1 + +) 


( — 9 ^ G * + + + ) sin (x - 2to)t 

\8x (x -u r 7 


3 kp ° r 

- — s — s-Tj — j cos (x+2io)t (1 4- 4- ) 

16(0 (x4-(o) (x -w ) L 

' C k . 1 %.a - 1 " 1 " 1 -) sin(x+2u>)tj 

8x (x -(0 ) 7 J 


+ 


9 kG^ 


; 2 2,4 

8 (x -to ) 


cos xt (1 4" 4-) -4 f 4-4^ sin xt 

8x(x‘ -w ) 7 .. - 


3 kG 1 


4 (3x 4- (j) (x4-(o) (x 2 -W 2 ) 3 


COS (2x4-(o)t (14- 4-) 
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a very 


-2 f 9kC j, t 2 - g + +) sin (2x+u)t 
ax (x -to ) 


3kG l 


2 2 3 

4 (3x-to) (x-to) (x -to ) 


cos (2x-io)t ( 1 + + ) 


c. 


9kG 2 t 


Q , 2 2.2 

8x (x -w ) 


+ +) sin (2x-u)t 


k G' 


A . 2 2,3 , 2 Q 2. . 

4 (x -to ) (x -9co ) (_ 


cos 


xt (1 + + ) -( 15 * +) smut 

^8x (x -to ) 


k G‘ 


„ 2,2 2.3 , 

32x (x -to ) L 


cos xt 


(1 ++) 


87 k G t 

. ,2 2.2 

•4x (x -to ) 


+ +) sin xtj 


§JL^_ — cos xt 

,.2 2,3 /Q 2 2. 

4 (x -to ) (9x -to ) L 


With the substitutions: 

.2 


9 k G 


q / 2 2.2 

8x (x -to ) 


= P 


(1 + + j sin **] 

v 8x (x -to ) 


(97) 


(98) 


G 


2 2 , 

x - to ) 


= a 


/ 2 2 , 

(x - to ) = m 


rough approximation for the first two terms gives 
f(t) = z |<cos wt - cos (x+p)tj 

. 3 T 9 . . cos 3 tot , 1 _ 

- k Z 5-—- COS tot + 9 9 + 2 

L 4m 4 (x 2 -9to 2 ) 32x 2 


(99) 

( 100 ) 


cos (3x+3(3)t 
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• ra n ;-!;) cos a-pw-p)* + tWuko 008 (x+2 “ +p,t 


r^r^T cos < 2x -“ + P )t 


4 (x 2 - On 2 ) 


/ , 15kz , , 
cos <x + -g^ ) t 

2 


1 ^ i , 87 kz 

- 2 — 003 < x + “45 > t 


32x. 


4{9x 2 -U 2 ) 


cos (x + 


■ 21 kz 2 
8x - 



( 101 ) 


Actually the number of terms available from computer results is 
still insufficient to establish Equation (101 ) as the approximate solution 
through the first two terms, (k° and k^). The assumptions in regard to 
the trigonometric relationships are obvious. Evidently one must choose 
suitable values of to, x, G and k. 

Figure 11 shows the usual numerical solution compared with the 
evaluation of Equation (101) for the parameters indicated. 



f(t) 



t seconds 


cn 

-a 
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11. DIRECT SERIES SOLUTIONS 


Direct series solutipns may be divided into two categories: 

Io Cases when f(t) is not only bounded but meets the requirements 
of Section 3 for the Laplace transform type of solution. 

Section 4 shows an example of the direct series solution for 
Equation (13) and Figure 1 shows the graphical result for tfie numerical 
compared with the analytic result for given numerical constants. 

II. Cases where f(t) may or may not be bounded. The transform 
method of solution does not apply but the derivative computer programs 
are identical for either Case I or Case II solutions. Case II solutions 
include all equations under Case I. 

The ,f open-end integral 11 of Section 3, Equations (4) and (18) may 
be interpreted as a defined operation which does not involve the Laplace 
transform. It does however provide a formal solution in the form of a 
Maclaurin's series. 

Define 


F(p) = p 



q 



n=0 


f n (0) 

n 


P 


( 102 ) 


The integration is to be carried out by parts and is, of course, a 
series of differentiations, q is the number of derivatives involved, 
usually less than 60. f n (t) is the n th derivative of f(t). 

Now if -i— is defined as in Heaviside 1 s calculus 


P 


1 


n 

P 


,n 




n ' 


(103) 


and Equation (102) 

«« - ZD 


is evaluated as 

f n (0)t n 
n ! 


n=0 


( 104 ) 



59 


15 

which is Maclaurin's series for f(t) through the q th term. Kaplan has 
examined equations of the type considered here and has shown that they 
always converge to an analytic solution for q— >>o. 

Figures 12 and 13 show results for summation and direct series 
solutions compared with the solution by numerical integration for two 
equations having Case I solutions. For the same number of derivatives 
the summation solution maintains accuracy for much larger values of time. 

Figures 14 and 15 show comparative results for problems under 
Case II. No summation solution is available in these two cases. 

Kaplan's method of showing convergence also allows an estimate 
for the remainder of the series after n terms so the accuracy can be 
estimated in solutions under Case II. 




. 60 



1.0 





f /, + (x+y)f / +xyf=kf 2 
k= -0.1, x = J . y=J_ 

/j sr 

f (o)^o, -r' (o)-i 

* Numerical Solution 

o o Analytic Solution -3 Terms 
x x Direct Series Solution 




















FIGURE! 13 






t seconds 








600 


500 
400 
300 
ZOO 
I 00 





3 


f"+(x+y)f / -t-xyf = kf 2 
k=Oi,x = - 1 , y ~ 1. 

7T /IT ' 

f(o)-0, f'(o) - I 
Numerical Solution 


x x Direct Series Solution 


figure: is 


4 5 

■*- t seconds 


7 


8 


64 


12J COMPUTER PROGRAMMING AND OPERATIONS 


While some simple non-linear problems may be solved by hand 
calculations, even relatively simple problems are best solved by the use 
of a digital computer. 

The programming is carried out using the same operations as 
would be used by hand calculations. Active memory is at a premium so 
the various operations should be broken down into what might be called 
irreducible steps. Results are stored on cards, paper tape or preferably 
magnetic tape. 

The numerical coefficients for the terms generated by the process 
of differentiation cannot be carried in floating point since the last digit 
of the coefficient is required. In fixed point one must carry coefficients 
involving some 112 bits or more and sign. A variable word length machine 
is desirable but not mandatory. , 

The format for the individual terms should be arranged to carry 
as many variables as practical in addition to the large coefficient since 
another problem may require additional variables. 

The basic programs required include a program for differentiation, 
one for summation, one for the numerical evaluation of the transforms and 
a program for the final evaluation of the time functions. 

The program for differentiation involves the input of the individual 
terms, the differentiation of this term to obtain the next term or terms 
and a compare and pack routine to store the generated terms in a list in 
memory. After all terms of a given derivative have been processed the 
next set of terms forming the next derivative is output to tape and input 
again to obtain the next derivative. Since each derivative depends upon 
the tape made from the previous derivative it is essential to check each 
operation carefully for errors if possible. 
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The summation operation can be performed using all available de- 
rivatives which will go into the computer at one time. 

The storage space required is considerably reduced if the de- 
rivative terms can be sorted on the basis of the power of the non-linear 
coefficient involved since this decreases the number of terms to be 
summed in one run. 

If the space required is still more than the available space in 
memory, additional sorting may be carried out on the powers of the initial 
conditions. 

The actual 11 summation" computer program is a multiplication 
operation. If the operational expression 

p + A 

.(p+B) (p+C) <p+D) 

is expanded in memory as a binomial expansion in negative powers of p 
and the whole list of terms is multiplied successively by (p+B), then (p+C) 
and finally by (p+D) the final result in memory is (p+A), the numerator. 

The numerators of the individual terms are then processed to 
include the number of roots in each operational expression. This will be 
called a " numerator tape " . 

A complete set of the required roots and the initial conditions plus 
the parameters are then input to an evaluation routine. The numerator 
tape is input, a term at a time, and for each term the machine evaluates 
the inverse Laplace Transform and stores the result with a compare and 
pack routine. The result is the "analytic solution" for the first few terms. 

While some operational solutions may be evaluated by hand, and the 
resulting time function summed into a closed form solution, the usual pro- 
cedure is to obtain the result as a numerical solution in the independent 
variable, time. This is then compared with the equivalent numerical so- 
lution to determine the range for which the analytic solution is accurate. 

The Heaviside shift program for the expansion type solution always 
holds data in the form , 
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(Numerator Polynomial) t n e 
(Denominator Polynomial) 

Except for the numerical coefficients all terms are in algebraic form. 
This program evaluates the solution faster than the derivative type of 
solution. Unfortunately the final result must be programmed to evaluate 
the function in terms of time. 

The direct series solution can be evaluated by substituting the 
numerical value of each constant into the series. Each constant is then 
multiplied by t n and added to obtain the sum of the finite series. 



67 


13. CONCLUSIONS 


The principal objective of this study was to develop suitable 
mathematical and programming techniques to obtain analytic solutions 
for appropriate types, of non-linear differential equations using a digital 
computer. 

Two general methods have been developed for the solution of the 
type of equation considered. For most bounded functions the expansion 
method is quicker and consequently less expensive in terms of computer 
time. The derivative method, however, is much more flexible since the 
solutions may be of bounded or unbounded functions. Summation solutions 
may be obtained, perhaps in several forms. (Under certain circumstances 
the results might, for instance, be summed as a Fourier Series.) While, 
the direct series type of derivative solution provides a numerical solution, 
it is very expensive in terms of computer time and is much less effective 
than the equivalent summation solution using the same data. 

The technique of checking the validity of equations through the use 
of numerical integration is a very important step in any analytic solution 
since it provides an engineering type of assurance that the " small para- 
meter, " k, is "small enough, n at least for that case, and that a sufficient 
number of terms for convergence is available. Figure 16 shows a case 
where a very large number of terms would be necessary for reasonable 
accuracy, It is surprising that the analytic and numerical solutions match 
as closely as they usually do since both are approximations. The numeri- 
cal solution has another important function, it points out any peculiarity in 
the solution which might otherwise have been overlooked. 

The appearance of secular terms in non-linear solutions has always 
been a handicap in hand solutions. Since these terms rarely seem to have 
any pattern it is pure drudgery to compute them by hand and accuracy is 
imperative. With computer solutions accuracy must be checked at every 
step, particularly in input and output, but over-all accuracy is more or 
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less assured. One may obtain a sufficient number of secular terms to 
approximate the function or one may be able to recognize the pattern and 
assimilate the terms into available functions. In either case the computer 
is an invaluable aid. 

This study may be regarded as of an experimental character. Much 
is yet to be done. No attempt has been made to make a thorough study of 
any of the equations considered. In many instances the same equation 
may have many different forms depending upon the value of the parameters 
and the initial conditions. Since initial conditions frequently appear as a 
part of the argument of the functions involved, the resulting solution may 
be quite complicated in structure. 

One may well ask why analytic solutions are of interest since 
digital computers may easily be programmed to provide numerical solu- 
tions. The answer is, of course, that numerical solutions can no more 
take the place of analytic solutions for non-linear work than they could 
for linear work. One analytic solution is likely to provide more under- 
standing of a given phenomenon than many numerical solutions. 

Many books and papers have been written on the theory of non- 
linear differential equations, but the student finds discour agingly few 
actual solutions. While the techniques which have been developed during 
this study are merely an entering wedge and appropriate for a particular 
type of equation, they point to the possibility of using a computer for 
analytic solutions in other non-linear areas and indicate that any real 
future progress in non-linear mathematics is likely to be closely 
associated with computer work. 
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- Ua3 0 ii e~yt(e“^ x ' } ' 2 ^ t -1) + (e^^H-l) 

( x-3y) (x+2y ) (x-2y ) ( 2x-y ) (x-3y ) (x-2y) 

- 2 a 3 c ke-y t (e- 2 xt -l) + 

xy(x-2y) y(2x-y) (x-2y) 

+■ ^c3 e -yt(e(y- 2 x)t, l) _ ^ a U c 3 e ~yt( e >(x+y)t, 1 ) 

(x+y)(2x-y)2 77 (2x-y)(x-2y)(x+y) 

t WtcV^Cefr- 2 ^-!) - 8a3o^e"y t (e-^ 2x4 '^ t -l) 

(2x.y)2(x-2y) xy( 2x4y) 

- - 8a3 c ^e“y t (e-fr^H-l) 

xy ( 2x-y) x (x-2y ) (x+y) 

+ 8a3 c U e-7 t (e (y ~ 2y)t -l) + l^cV^ (e“( x+2y)t -l) 
x(x-2y)(2x-y) y(x-3y) (x+2y) 

- - 8a 2 c?e“yt(e“^^-l) ‘ 

y(2x-y)Xx-3y) y(x+y} (x-2y) 

+ Uag^e-^Ce^^^-l) - h^c 2 ^ (Jy-^H-l) 
y(2x-y) 6c-2y) x(Ux-y)(3x-y) 

+ ^Ve-^Ce^-l) + Ua^^e-ytCe^y-^H-i) - i+a^cVy* (e^-l) 
x 2 ( 3x-y ) x(3x-y)(2x-y) x 2 (2x-y) 


Cont. 
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B0 


- (e^ 3 *^-!) + Ua^c 3 e ( e -xt -1 ) + 8a^c 3 e~ yt; («^ y ~ 3x ^ t '-l) 

~~ 3ac 2 y x 2 y y(3x-y)(2x-y) 

- Ma^e^e^-l) + Ua3c^e“y t (e-^ 2x4 y )fc -l) - 1 

xy( 2x-y) (x-2y ) ( 2x+y) (x+y) ' x (x-2y ) (x+y) 

- + 8 a^c^e“^'(e” x *’-l) 

(3x-y ) ( 2x-y ) (x-2y ) x(2x-y) (x-2y) 

- iia^c^e-y^ ( e ^ -1 ) + iia^c3e~^(e“ x k.-l) + lia^c^e**^ (e*^^- 1) 

(iix-y)(3x-y)(2x-y) . x(3x-y)(2x-y) x^(2x-y) 

- U ti c 3 e“y t (e~ xt -.l) - ha^cK^ (e" 3 ^-!) + l^cVy^e" 3 ^-!) 

-2^ 2x _-)- ; 3x3 

+ ^c^e-^Ce^-l) - 8a 3 c^e" , y t (e~ :xt ~l) - ^cV^Ce^ 2 *^-!; 
x 3 x 3 y(2x+y)(x-*y) 

+ ^e^e-yfcte-^-l) + 2a 2 A ^(e- 2 ^-!) - ^aVe^e^-l) 
xy(x+y) x 2 y x 2 y 


- 2a^c 3 e*y t (e“ 3xfc -l) + 2a U cVy t (e^-l) + (e” 2xfc -l) 


- kaVe-ytCe-xt-l) - 8a3cVy t (e-( 2x4 y) t -l) + 8a3 c Vyt(e‘ xfc -.l) 

x3 y(2x+y)(x+y) xy(x+y) 


+ i;a 3 c^e*“^ (e** 2x ^-l) - l6a 3 c^e~y^ + 2a 2 c^e**^ 

x 2 y x 2 y y(x-2y) (x+2y) 

+ - 2 a 2 c g e“ - y t (e~' 2 xt ~l) + (e^-l) 

xy(x-2y) x^(x-2y) x 2 (x-2y) 

- 2a 3 c^e~yt(e“ 3xt -l) + 2a 3 c^e-y t (e~ xfc -l) + Ua.^e'"^ . 

__ 3x 2 (2x~y) x 2 (2x-y) y( 2x-y)(x4y) 
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- Iia^c^e"^ (e" 3 ^-!) - 8a^o^e“yk(e"(^ x ~ty}t-l) + Sa^c^e"*^ (e** 3 ^-!) 

xy(2x-y) x(2x+y) (x+y) (x+y) 

H - aacV^g-fr^yH-l 

xy(x+y) y 2 (x+2y) " 

+ ZacV^Ce-^-l) + Uac^e~y^(e~^ x+ y^~l) - UacV^Ce-**-!) 

xy 2 y 2 (x+y) xy 2 - 

- 2aVe-y t (e-3xt_ 1) + ga^oV^e^-I) + 2a5c 2 e-yt(e- 2xt ..l) 

3x3 x 3 

- UaVe-^Ce-^-l) - + 8a^ C 3e~yt(e~ xt ~l) 

x 3 y(2x+y) (x+y) xy(x+y) 

- ^cV^fe^-l) + 2a3 c ^e“y t ( e < -^ x+2 y) t -l) - 2a3c U e^ t (e^ fc -l) 

~ ?y y(x-2y)(x+2y) xy(x-2y) 

- 2a3 c Vy t (e- 2xt -l) + (e^-l) - ga^cV^Ce'^-l) 

x 2 (x-2y) x 2 (x-2y ) 3x 2 {2x-y) "3 

+ 2a i; c3e"y fc (e- xfc -.l) + UaMe-y^e-fr*^-!) 

x^(2x-y) y(2x-y) (x+y) 

- 8 a 3 c ^e“-y t (e- > ^ 2x4 y) t -l) + 8a3A“ yt (e^-l) 

x(2x+y)(x+y) x 2 (x+y) 

+ 8 a 3 c i+ e“y t (e-^) t -l) - 2a 2 c g e"y t (e~ (y+2 y^ t -l) 

y 2 (x+2y) 

+ + l^A^Ce^^l) - l»2oS»-y t (e'* b -l) 

xy 5 y^(x+y) xy 2 

- Ua^c^e~^ rfc '(e"'^ 2x '^ r ^~l) + Ua^c^e**y b (e* >:ict> -1) 

- x(2x+y)(x+y) x^(x+y ) ' 


Cont. 



xy(x-ty) y 2 (x+2y) 


+ + 8a3 c ^ e -yt(e-^ + y) t ~l) - 8a3cV^£^0 

xy^ y2(x+y) xy 2 

- Ua 2 c^e"yt(e-3yt^ 1 )+ Ua 2 c^e“y t (e- Xt -1) + ka 2 c£ e * - J rfc (e“ xfc -l) 
_ 3y(x-3y) (x-2y) x(x-3yKx-2y) ~ xy(x-2y) 

- ^c^e-^Ce-C 2 ^)^-!) + Ua3 c ^ e -yt( e -xt_ l) . 2a3 c Vy t (e~ 2 yt..l) 

(2x+y) (2x-y) (x+y) x(2x-y)(x+y) y(2x-y)(x-2y) 

- + h3 2 c S e-V t (e- :x *'-l) - ^cV^Ce- 2 ^-!) 

xy(x+2y) x 2 y xy(x-2y) 

+ ^ac 6 e-y t (e~3yt-i) - hac^e^ (e^-l) - 2ao 6 e^ t (e^-1) 

3y2 (x-3y) xy(x-3y) y2(x-2 y) 

+ Uac^e~y^ (e - ^-!) - l4a^c^e“ , y* : '(e"’3 :iC ^_i) + l;akc3e~y^(e”^--l) 
xy(x-2y) . 3x' 2 (3x-y) xy(3x-y) 

- + ItaVe^Ce^-l) + ( e “ 2xt -l ) 

xy(2x-y) xy 2 xy(2x-y) 

- 8a3 c ^e-y t (e-y fc -l) + ^cV^e-fr 42 ^-!) - Ua^e-y^e"^-!) 

y 2 (2x-y) (x+2y) (x-2y) (x+y) y(x-2y)(x+y) 

- 2a 2 c g e“y t (e~ 2xt ~l) + (e-yt-1) - Ua^e^Cg^f^ 

x(2x-y) (x-2y) y (2x-y)(x-2y) 3x(3x-y) (2x-y) 

+ Ua3 c k e --yt(e-yt.-l) + Ua3 c ^e~y b (e“ (x+y)t -l) - 

y(3x-y)(2x-y) x(2x-y) (x+y) xy(2x-y) 

- lia 2 oge-y t (e-^ 2x4 y) t -l) + Ua^e^e"^) 

x 2 (2x+y) x 2 y 
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- ^ a6 ce**yfc ( e (y-*^PcH»i) + 2a 6 ce"y*(e(y-2x)t-i) 

x 2 (ipc-y) x^(2x-yj 

+ ij.a^ce“y^(e^y“3x^-l) — lj.a^ce”’ y k.(e( y “ 2x ^— l) 
x 2 (3x-y) x 2 (2x-y) 

- 8a^’c 2 e~ y k (e"^ 3 ^"!) + 8a 3 c 2 e~ y k(e^ y ~ 2x ^-l) 

3xy (x+y ) y(2x-y) (x+y) 

+ Sa^cV^Ce^- 3 *^-!) - lgaVe^Ce^^-l) 

xy(3x-y)' xy(2x-y) 

+ 2 a K%^ (e“ (2x4y)t ~l) - 2aU c 3 e ~yt 

y(2xty)(x-2y) y(2x-y)(x-2y) 

- W^V^Ce^" 3 *^-!) + jAV^ (e ( y - 2 *H..l) 

xC3x-y)(x-2y) x(2x-y) (x-2y) 

- 2aVe~ yfc (e( y ^ x ) t ~l) + 2a5c 2 e" yt (e( y - 2x > t -l) 

x(ipc-y) (2x-y) x(2x-y) 2 

+ 2a^c 2 e~y t (e- 2:)ct -l) - UaVe"^ (e( y - 2x H~l) - 8a^c 3 e“y t (e~3xt^ 1 ) 
xy(2x-y) y(2x-y)2 3x 2 (x+y) 

+ 8a^c 3 e"’ y k (e (y _2x )^-l) + 8a^c 3 e -y ^(e“ 2x ^-l) 
x(2x-y) (x+y) x 2 y 

- 8a^c 3 e~‘ y ^(e^ y ~‘ 2x ^ t ~l) - 2a.^e^ ( e -( 2 *+ y )t-i) 

xy(2x-y) y 2 (2x+y) 

+ 2a 3 oV yfc (e( y ~ 2x > t -l) + 2a3 c ^- yt (e^-l) 
y 2 (2x-y) xy 2 

- Ua3 c ^e~ yfc (e( y “ 2x ) t -i) - UaVe"^ (e" 3 ^-!) 

y 2 (2x-y) 3x 2 (x+y) 


Cont. 
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8.4 

T T8a 2 c^ e -y t ( e -^ )t -l) - 8a 2 c^e"y t (e“y t ~l) 

x 2 (x+y) x^y 

- Itac^e**^^ (e’( x+g y ^-1) + liao^e*"^ (e~y^~l) 

y(xty)(x+2y) y 2 (x+y) 

+ kac^e~^(e~^ y4y ^‘-l) - 12ac^e'~ y ^ i (e~^ rb ~l) 
xy(x+y) xy 2 ™ 

- 2a3 c 1 ^e-yb( e -(2x4y)t_ I ) + 2a 3c^ e -yt( e «-yt_i) 

x 2 ( 2x+y ) x 2 y 

“ - lia3A-yt (e^-l) 

x 2 (x+y) x^y 

- 8a 2 c^e“’3 r ^ (e“( x+2 5 r ^-1) + 8a 2 c^e“^(e**^— 1) — 8a 2 C' l ’e“‘’^'(e“’^ , -l) 

x(x+y)(x+2y) y 2 (x+y) xy 2 

+ gac^-ytCe^yt^x) _ 2 ac 6 a -yt( e -.yt.x) . UaA"^ (e^* 4 ^-!) 
3y^"(x--2y) y^'(x-2y) x(x +y)fc- 2y) 

+ itac^e -yt ( e "3^-1 ) - 2a 2 c^e" yt (e~^ 2x4y ^-l) + 2a 2 c5e~y t 
xy(x-2y) x(2x+y)(2x-y) xy(2x-y) 

+ 2a 2 c 5 e-yt(e" 2 y t ~l) - Ua^e^Ce^-l) - BacV^e-fr^*-!) 

y 2 (2x-y) ~ y 2 (2x-y) x(x+y)(x+2y) 

+ SaA^Ce-yt-l) + liaA-ytCe^yfc-l) - 2cVyt(e"3yt,i) 
xy(x+y) xy 2 3y3 

+ 2c ? e“y t (e-yt-i) + 2cVy t (e- 2 y t -i) - kcV^Ce-yfr-i) 


plus terms having a multiplier with a power of k larger than 3. 
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